home *** CD-ROM | disk | FTP | other *** search
/ OpenGL Superbible (2nd Edition) / OpenGL SuperBible e2.iso / tools / GLUT-3.7 / LIB / GLE / UROTATE.C < prev    next >
Encoding:
C/C++ Source or Header  |  1998-08-12  |  6.2 KB  |  218 lines
  1.  
  2. /* 
  3.  * MODULE: urotate.c
  4.  *
  5.  * FUNCTION:
  6.  * This module contains three different routines that compute rotation
  7.  * matricies and return these to user.
  8.  * Detailed description is provided below.
  9.  *
  10.  * HISTORY:
  11.  * Developed & written, Linas Vepstas, Septmeber 1991
  12.  * double precision port, March 1993
  13.  *
  14.  * DETAILED DESCRIPTION:
  15.  * This module contains three routines:
  16.  * --------------------------------------------------------------------
  17.  *
  18.  * void urot_about_axis (float m[4][4],      --- returned
  19.  *                       float angle,        --- input 
  20.  *                       float axis[3])      --- input
  21.  * Computes a rotation matrix.
  22.  * The rotation is around the the direction specified by the argument
  23.  * argument axis[3].  User may specify vector which is not of unit
  24.  * length.  The angle of rotation is specified in degrees, and is in the
  25.  * right-handed direction.
  26.  *
  27.  * void rot_about_axis (float angle,        --- input 
  28.  *                      float axis[3])      --- input
  29.  * Same as above routine, except that the matrix is multiplied into the
  30.  * GL matrix stack.
  31.  *
  32.  * --------------------------------------------------------------------
  33.  *
  34.  * void urot_axis (float m[4][4],      --- returned
  35.  *                 float omega,        --- input
  36.  *                 float axis[3])      --- input
  37.  * Same as urot_about_axis(), but angle specified in radians.
  38.  * It is assumed that the argument axis[3] is a vector of unit length.
  39.  * If it is not of unit length, the returned matrix will not be correct.
  40.  *
  41.  * void rot_axis (float omega,        --- input
  42.  *                float axis[3])      --- input
  43.  * Same as above routine, except that the matrix is multiplied into the
  44.  * GL matrix stack.
  45.  *
  46.  * --------------------------------------------------------------------
  47.  *
  48.  * void urot_omega (float m[4][4],       --- returned
  49.  *                  float omega[3])      --- input
  50.  * same as urot_axis(), but the angle is taken as the length of the
  51.  * vector omega[3]
  52.  *
  53.  * void rot_omega (float omega[3])      --- input
  54.  * Same as above routine, except that the matrix is multiplied into the
  55.  * GL matrix stack.
  56.  *
  57.  * --------------------------------------------------------------------
  58.  */
  59.  
  60. #include <math.h>
  61. #include "gutil.h"
  62.  
  63. /* Some <math.h> files do not define M_PI... */
  64. #ifndef M_PI
  65. #define M_PI 3.14159265358979323846
  66. #endif
  67.    
  68. /* ========================================================== */
  69.  
  70. #ifdef __GUTIL_DOUBLE
  71. void urot_axis_d (double m[4][4],         /* returned */
  72.                   double omega,         /* input */
  73.                   double axis[3])        /* input */
  74. #else
  75. void urot_axis_f (float m[4][4],         /* returned */
  76.                   float omega,         /* input */
  77.                   float axis[3])        /* input */
  78. #endif
  79. {
  80.    double s, c, ssq, csq, cts;
  81.    double tmp;
  82.  
  83.    /*
  84.     * The formula coded up below can be derived by using the 
  85.     * homomorphism between SU(2) and O(3), namely, that the
  86.     * 3x3 rotation matrix R is given by
  87.     *      t.R.v = S(-1) t.v S
  88.     * where
  89.     * t are the Pauli matrices (similar to Quaternions, easier to use)
  90.     * v is an arbitrary 3-vector
  91.     * and S is a 2x2 hermitian matrix:
  92.     *     S = exp ( i omega t.axis / 2 )
  93.     * 
  94.     * (Also, remember that computer graphics uses the transpose of R).
  95.     * 
  96.     * The Pauli matrices are:
  97.     * 
  98.     *  tx = (0 1)          ty = (0 -i)            tz = (1  0)
  99.     *       (1 0)               (i  0)                 (0 -1)
  100.     *
  101.     * Note that no error checking is done -- if the axis vector is not 
  102.     * of unit length, you'll get strange results.
  103.     */
  104.  
  105.    tmp = (double) omega / 2.0;
  106.    s = sin (tmp);
  107.    c = cos (tmp);
  108.  
  109.    ssq = s*s;
  110.    csq = c*c;
  111.  
  112.    m[0][0] = m[1][1] = m[2][2] = csq - ssq;
  113.  
  114.    ssq *= 2.0;
  115.  
  116.    /* on-diagonal entries */
  117.    m[0][0] += ssq * axis[0]*axis[0];
  118.    m[1][1] += ssq * axis[1]*axis[1];
  119.    m[2][2] += ssq * axis[2]*axis[2];
  120.  
  121.    /* off-diagonal entries */
  122.    m[0][1] = m[1][0] = axis[0] * axis[1] * ssq;
  123.    m[1][2] = m[2][1] = axis[1] * axis[2] * ssq;
  124.    m[2][0] = m[0][2] = axis[2] * axis[0] * ssq;
  125.  
  126.    cts = 2.0 * c * s;
  127.  
  128.    tmp = cts * axis[2];
  129.    m[0][1] += tmp;
  130.    m[1][0] -= tmp;
  131.  
  132.    tmp = cts * axis[0];
  133.    m[1][2] += tmp;
  134.    m[2][1] -= tmp;
  135.  
  136.    tmp = cts * axis[1];
  137.    m[2][0] += tmp;
  138.    m[0][2] -= tmp;
  139.  
  140.    /* homogeneous entries */
  141.    m[0][3] = m[1][3] = m[2][3] = m[3][2] = m[3][1] = m[3][0] = 0.0;
  142.    m[3][3] = 1.0;
  143.  
  144.  
  145. }
  146.  
  147. /* ========================================================== */
  148.  
  149. #ifdef __GUTIL_DOUBLE
  150. void urot_about_axis_d (double m[4][4],         /* returned */
  151.                         double angle,         /* input */
  152.                         double axis[3])        /* input */
  153. #else
  154. void urot_about_axis_f (float m[4][4],         /* returned */
  155.                         float angle,         /* input */
  156.                         float axis[3])        /* input */
  157. #endif
  158. {
  159.    gutDouble len, ax[3];
  160.  
  161.    angle *= M_PI/180.0;        /* convert to radians */
  162.  
  163.    /* renormalize axis vector, if needed */
  164.    len = axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2];
  165.  
  166.    /* we can save some machine instructions by normalizing only
  167.     * if needed.  The compiler should be able to schedule in the 
  168.     * if test "for free". */
  169.    if (len != 1.0) {
  170.       len = (gutDouble) (1.0 / sqrt ((double) len));
  171.       ax[0] = axis[0] * len;
  172.       ax[1] = axis[1] * len;
  173.       ax[2] = axis[2] * len;
  174. #ifdef __GUTIL_DOUBLE
  175.       urot_axis_d (m, angle, ax);
  176. #else
  177.       urot_axis_f (m, angle, ax);
  178. #endif
  179.    } else {
  180. #ifdef __GUTIL_DOUBLE
  181.       urot_axis_d (m, angle, axis);
  182. #else
  183.       urot_axis_f (m, angle, axis);
  184. #endif
  185.    }
  186. }
  187.  
  188. /* ========================================================== */
  189.  
  190. #ifdef __GUTIL_DOUBLE
  191. void urot_omega_d (double m[4][4],     /* returned */
  192.                    double axis[3])        /* input */
  193. #else
  194. void urot_omega_f (float m[4][4],     /* returned */
  195.                    float axis[3])        /* input */
  196. #endif
  197. {
  198.    gutDouble len, ax[3];
  199.  
  200.    /* normalize axis vector */
  201.    len = axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2];
  202.  
  203.    len = (gutDouble) (1.0 / sqrt ((double) len));
  204.    ax[0] = axis[0] * len;
  205.    ax[1] = axis[1] * len;
  206.    ax[2] = axis[2] * len;
  207.  
  208.    /* the amount of rotation is equal to the length, in radians */
  209. #ifdef __GUTIL_DOUBLE
  210.    urot_axis_d (m, len, ax);
  211. #else
  212.    urot_axis_f (m, len, ax);
  213. #endif
  214.  
  215. }
  216.  
  217. /* ======================= END OF FILE ========================== */
  218.